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ABSTRACT 


Supersonic flow past harmonically vibrating two- 
dimensional panels and cylindrical shells is analyzed 
using a linearized method-of-characteristics procedure. 
A detailed description of this method to solve the 
linearized unsteady potential equation is given and 
numerical results are presented to indicate the nature 
of the aerodynamic pressure ds а Also, com- 
parisons of the present results are made with earlier 
work by Nelson and Cunningham and by Anderson which is 


based upon quite different approaches. 
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TI. INTRODUCTION 


The determination of the pressure distribution over | 
vibrating panels and shells is a prerequisite for the E 
of aeroelastic stability. This problem has been previously 
considered by a number of investigators, who applied either 
various asymptotic approximation techniques (e.g., piston 
theory) or operational methods to the solution of the lin- 
earized unsteady potential equation. For a recent compre- 
hensive review of this subject refer to Reference 7. In 
the present study a quite different approach is developed; 
1.е., the method of characteristics is used to obtain a 
solution to the problem of supersonic flow past vibrating 
panels and shells. 

The major assumptions in the mathematical model are 
that we are dealing with a perfect gas, that the flow is 
irrotational and supersonic, and that linearization is 
permissible. Reference 1 was used to help develop tne 
linearized unsteady potential equation and the boundary 
conditions for panels and shells. These equations were 
then rewritten in cylindrical coordinates. 

Steady supersonic flow past bodies of revolution was 
studied by several authors using linearized characteristics 
methods; i.e., Haack [Re£. 9|, Sauer [Ref. 10], Oswatitsch 
and Erdmann [Re£ . 11. Teipel [Re£. 12] developed a char- 
acteristics procedure for two-dimensional supersonic flow 


past airfoils oscillating at arbitrary frequency using the 





velocity perturbations and the velocity of sound perturbation 
as dependent variables. Platzer and Sherer (Ref. 13] gave an 
extension of the Oswatitsch-Erdmann method to slowly oscil- 
lating bodies of revolution. 

In the present report a linearized characteristics method 
is presented to study supersonic flow past cylindrical shells 
vibrating at arbitrary frequency. The velocity potential 
function is introduced as the dependent variable and the 
linearized unsteady potential equation is solved subject to 
the proper boundary conditions. Although two-dimensional 
flow past vibrating panels is contained in this solution as 
a Special case for infinite shell radius, two separate 
programs were written in Fortran IV code for the IBM 360 
манс, “опе for panels and one for sielils. fits, слем NEM 
program is complementary to Teipel's work in that it uses 
the velocity potential rather than the velocity perturbations 
as dependent variables. The method of singularities [Re£. a] 
served as an independent method to verify the results. 
Anderson's work, on the other hand, was used to compare the 


method-of-characteristics predictions for cylindrical shells. 








II. PROBLEM FORMULATION 


In nearly all engineering studies some basic assumptions 
have to be made to formulate the problem. In this study l 
this assumption 15 that the compressible medium with which 
"ме are dealing is a perfect gas. We therefore have the 


equation of state for a gas 


P = pRT (ШІ) 

Flow of a perfect gas is then completely described bç 
expressing the ambient quantities P,p and T as well as the 
velocity components u, v, w as a function of position x, y, 
z and time t. These six dependent variables are always 
related by the conditions of conservation of mass, momentum 
and energy as developed in Johns [вее. Це 
Conservation of Hass 

nn en. Eon ОЕА 


Conservation of Momentum 


I 2 _ ач, 5 38 + vot + wd 

p 9х dt 9х ду 92 

і ЭР _ av, | av aV ЗУ 

р а с о Ws (2.3) 
A и оа 
EE ааа oz 

/ 
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Conservation of Energy 


D w? 1|9 ,— 2 3 
Dt E = => 3x (МР) zh ay (VP) т уж (Р) (2.4) 


where _ 
W2 = u2 + v2 + w?. (2.5) 


Equation 2.4 can also be stated as (See Ref. ]) 


ne 
DE 0 | (256) 


The condition of conservation of energy therefore reduces 

to the simple statement that the entropy of a fluid particle 
remains constant. For later convenience we now rewrite the 
continuity equation, Eq. 2.2 


Do да | ove ow) _ 
P+ + у + = о | Мана 


Introducing the speed of sound 


a? = 35 (2.8) 
we have 
oy fa me SL pe 
t ӘР DE > Dt 
Therefore the continuity equation can also be written 
= . (2 + 57 +32) (2.9) 


The velocity potential is introduced to reduce the 
number of dependent variables. Its existence depends upon 


the condition of irrotationality of the flow. This condition 
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is expressed mathematically by the vanishing of the curl 


of the velocity vector and can be written in component form 


а? 
< 
а? 
el 
Qo 
z 
а» 


У да ди г 
3x - 37 50, DONE. 0,52 = COE (2.10) 
It can be shown by Kelvin's theorem кес. 2 that flows 
starting elther from reservoirs or from parallel stream- 
lines and being irrotational initially will remain so. Only 
strong shock or intense heat will require a reformulation 

of this law. However, such effects are usually quite 
unimportant in aeroelastic considerations. Therefore, Eq. 
2.10 is necessary and sufficient to derive the velocity 


components u, v, w from the velocity potential $ such that 


ч = 2% y E = ‚м E | (2.11) 


Therefore the continuity equation can be written 








| 
т) 
о 
a 
с 
а? 
< 
а? 
2, 
а? 
~ 
X 
а» 
~ 
+> 
+ 
Qo 
N 
ме 


(2.12) 


о? 
„а 
о? 

< 
@ 
N 


In vector form the Eulerian equations of motion can also 


be written as 


-> 


QV E u Т 
асва = = 3 grad P. (220595 


This is equivalent in irrotational barotropic flow to 


д w? Т | 
sg grad y + grad (xA w grad P (2.14) 
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Or 
oo wt, [а 
grad SF oe T D = 0 


Integration therefore gives ES 


оф 
de + г. |ж. C (t) (2515) 


Moreover, C(t) - l u? whenever the remote fluid motion 
2 
consists of parallel streamlines with ale U. The 


partial derivative of Eq. 2.15 with respect to time is 


д эф мА 3 АР _ 1 ӘР 
- 3€ 2] - = р о 38 — 


In rewriting the continuity equation 


Я M ES ЈЕ 3P ЭР, ЭР] _ 
pes utt —— С 2 5535452 й о 20 
а2р 


one obtains therefore as a result of Eqs. 2.3 and 2.16 


$° e? 2 
1- |, +)1-—] $ а NE 
a2 XX ED уу m 22 
О, 246 204 
M X'y 2 Y 2 
Шы 'xy ^ uz іка” “а був" — 
2%. 20, 26, I 
“аз xt 7 22 a nos 


Let a denote the speed of sound. For infinitely large 
speeds of sound (incompressible flow) the above equation 


= 
r 


reduces to Laplace's equation. Equation 2.18 is the 
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governing differential equation for general nonstationary, 
nonviscous, potential flow and is not limited to small 
disturbances. However, because of its nonlinear nature, 
solutions were found only in very few special cases.  There- 
fore, it is necessary to introduce the small disturbance 
concept. Regarding all velocity disturbances small in 
comparison with U, a, and (U - a) and pressure differences 
and density changes small in comparison with main stream 
pressure and density, a perturbation potential Y is introduced 
so that 

ф = Ох + Фф (2.19) 
If a is expanded around its value for the undisturbed stream 
the procedure of retaining only first order terms leads to 


the following equation for the disturbance velocity potential 


(1-M?)f ... +P yy + ZZ - =? xt ^ 5 Фе ze) (2.20) 


where c is the constant value of the speed of sound in the 
uniform stream and M the constant free-stream Mach number. 

It is this linearized unsteady potential Eq. 2.20 which is 
being used in most aeroelastic analyses as the basic equation 
to predict the oscillatory pressures and een aero- 
dynamic forces. However, it must be emphasized that its 
validity is restricted to purely subsonic or supersonic 


flows, as well as to sufficiently unsteady transonic flows. 
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Expressing Eq. 2.20 in cylindrical coordinates one has 


«Ф I 1 Фф 1 2p _ 
(2.21) 


There are two main boundary conditions applicable to a 
vibrating cylinder in supersonic flow. The first is that 
the panel surface is impenetrable to the gas, i.e., the flow 
is tangential to the surface at each instant of time. The 
second condition is prescribed by the Sommerfeld radiation 
condition which requires waves to propagate away from sources 
of disturbance toward infinity. Furthermore, uniform super- 
sonic flow must exist upstream of the panel leading edge, 
since no disturbances can propagate upstream. 

If the equation of the surface of a body which is moving 


in a time-dependent manner is given by 


E(x Sy 2 í; t) a0 (2.22) 


then it can be shown that the flow tangency condition 


requires [Ref г J 


DF ЗЕ Е Е oF Е 
Dt ^ 35€ * V 3x *t V у += 0 (2454 


This condition states that the rate of change of the numerical 
value of the function F is zero, when we follow the motion 

of a particular fluid element (substantial derivative), SO 
that the element continually touches the surface F - 0. 


Consider a two-dimensional panel (Fig. la) whose chord is 
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aligned with the x-axis and whose upper surface is exposed 
to a uniform flow of velocity U in the positive x-direction. 


We then find that we can describe the body by 


= 


ee Ў (2.24) 
Since 3Е = 1 we have from the boundary condition 
92 
ah — oh _ 
A ае lone gn cue (2 226) 


Imposing the linear perturbation assumption Eq. 2.19 we have 


d = Ux +P (2.19) 


where the disturbance velocity components are obtained from 


na 9% мол оф (2.26) 


эх ' 92 


and satisfy the order-of-magnitude requirement 


Ж, W << U . (2.27) 


шые terms u = can then be neglected compared to the much 
oh 


larger term . This leads to 
. àh ah 
м = EU (2.28) 


It can be shown that it is possible to express the flow 
tangency condition in the plane z = 0. Taylor series 
expansion gives 


3w(x,0*,t) 


w(x,z,t) = w(z,0*,t) + h Z 


К + 
үш То ох О ЛЕЛ з. (2.29) 
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FIGURE 1b PROFILE AND CROSS SECIICN 


OF CYLINDRICAL SHELL 
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Using the same argument as above the product terms are 


neglected so that the flow tangency condition equals 


m 3h + | 
w > + U SE for z = 0 (2.30) 


A similar consideration for the cylindrical shell (Fig. lb) 


whose vibration mode is given by 


iot 
h(x,0,t? = Z(x) Cos n9 e (231) 


results in the following boundary conditions 


: lot 
Ar ax + іші Cos n9 e 


r-R 


g) о [зё 


for O < x < L 


= 0 (2.32) 
for 7x 250 


The calculation of perturbation pressures and pressure 


coefficients begins with Eq. 2.15 


9 
д 


<> 


W 2 dP 
a Т z 7 C (t) 


rr 


For uniform parallel flow far upstream we have 


172 
C(t) = 2 


Using again the perturbation assumption, Eq. 2.19 
ф = Ux +P 


and noting that 


2 2 2 
2 = af А эф 
W UI N уі ава а (2.5) 
We can approximate 
| му W жол US оке (2588) 
ot О die 2 9х 





by dropping second order terms. Also the density-pressure 


relation is expanded as 











В 2 5 ” 
p (p) = p + de = 7 і. а“р å 
со ар (P г.) 2! ар? р P TERCERA (2.34) 
Therefore 
p p p 
dp m dp _ “1 ' р-р. Р-Р. 
p (p) do а D 7 рү J89 З р. 
P a 1 +\ap ps 
Ро Bo m p 
(2,35) 


Hence, the perturbation pressure is related to the perturbation 


potential by 


д д 
Р-Р. = - DM E: + Џ S= (2.36) 


Writing this equation in pressure-coefficient form one has 


РР _ . 2. 90 _ 2 29 
= Gee ae jae (2 37) 


C 


Consistent with the assumed cylinder deflection, Eq. 2.3l a 
velocity potential is now assumed in the form 


lwt 


q = $(x,r) Cos n8 e (2.38) 
Substituting this into the Eq. 2.21 results in 
2 21090) 2 
Eo 1 n 10) ü) HE 
(1-М4) Ф sls Фү + = 2 E pa EX Фу T zo 0 (2.39) 


This is the basic equation for é(x,r), which needs to be 


solved, subject to the boundary conditions 


o 
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ге 24: 
P SUE 104 ога аа 
r=R 
= 6 For x < 0 | (2.40) 
| = Ж - T 
When non-dimensional coordinates X - r, r з г and the 
wL | ; - 

reduced frequency K = u are introduced then yt Tszreund 
(1-M2) 6 + 6 + 2 Фо - п? 6- 21KM26 + K2M26 = O (2.41) 

хх A x 4 


The bar has been omitted for simplicity. The boundary 


condrtion is 


қ 
|| 


Пра 
24 + iz O < x < L 


х < 0 (2.42) 
= 0 


Converting Cp to cylindrical coordinates and using the above 


non-dimensional coordinates one has 


Ср = -2 Е _ Я Cos née iKt (2.43) 
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III. NUMERICAL TECHNIQUE 


A. METHOD OF CHARACTERISTICS 


Consider the following pair of simultaneous first order 


partial differential equations for the dependent variables 


u and v, 


1 эх 97 J. EAE 

gu gu ду У 
А. о ЗОЯ No 
where the coefficients Ај, МИРА D. 


and v but not of x and y. 


0 (321) 


are functions of u 





ENITIRL Сом огол 


FIGURE 2 
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Assume that in Figure l the solution is known up to the 
curve CPC. At P one knows continuously differentiable 
values of u and v along the curve CPC as well as all their 
derivatives in the directions that lie below the curve. 
However, the question needs to be answered whether there is 
enough information to deduce the directional derivatives 

of u and v at P in directions that lie above the curve. 


The directional derivative Of u in the directions Tis: 


ou du ox du ду 


95 ox 95 Эду IS (3.2) 
so that the directional derivative in any direction is known 
as soon as the derivatives 3u/90x and du/dy are known. 

To find the values of 93u/3x and du/dy, write out 
equations 3.1 at point P and the directional differentials? 


of u and v taken along CPC at P. In matrix form the four 


equations are: 





Ај By Ci 9u/9X 0 
А. B> Co эа/ду 0 
ах als 0 IV/IX = Wu (3.3) 
- 0 0 ах ом/ду ES 
With u and v known at P the coefficients Ај, .. 77 р, 


will also be known. When the directions of CPC are known 


then dx and dy can be found. If u and v are then known along 


ki Eg. 3.3 du is ап abbreviation La 91/95 ds, where 
9u/9s is given by Figure 3.2 and s now measures distance , - 
along the curve CPC. The symbols dv, dx, and dy stand for 
similar abbreviations. 
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CPC, du and dv will be known. Equations 3.3 therefore 
constitute a pair of four simultaneous linear algebraic 
equations for the four first derivatives. If the determinant 
of this matrix is not zero, there is a unique solution. In 
this case satisfaction of these governing equations requires 
that the directional derivatives have the same value above 
and below the line CPC. If the matrix determinant had 
equaled zero, there would not be a unique solution and 
there could be discontinuities across the line CPC. 

Expanding the determinant of Eqs. 3.3 and setting it 


equal to zero one has 


= Pape мы = 
Es А.С1)4у (А.Б, A4D*B4C5 BoC, )dxdy + 


'(B1D2-B5D,)dx? = 0 | (3.4) 


This is a quadratic equation for the slope dy/dx. If the 
direction of the curve CPC at P is such that it has a slope 
satisfying Eq. 3.4, then the derivatives of u and v are not 
uniquely determined by the values of u and v along the curve. 
Such a direction is called a characteristic direction. The 
quadratic Eq. 3.4 gives two real slopes, one real slope or 

a pair of complex slopes, depending upon whether the 


висит папе 


(A,D4-A4Dj*B 


= аа y = 
192 В2С1) о А-С1) (В10- BoD, ) (3-5) 


is positive, zero or negative. This is also the criterion 
for cataloging the system Eq. 3.1 as hyperbolic, parabolic, . 


or elliptic. The system is hyperbolic at a point if there 
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are two real characteristic directions; it is parabolic if 
tere is only a single characteristic direction sorge 
elliptic if there are no real characteristic directions at 
the point. 

The foregoing analysis can be repeated for the single- 


order quasi-linear equation, 


д 2 oe 32 
а 5-7 + E С 55 = $ (3.6) 
е е l е QU 
ШИ ОС а, Б, c, and £t are functions of x, y, 3, ee ome 


Әш/9у. The characteristic directions are determined from 


the quadratic 


ady? - bdxdy + cdx? = 0 (3.7) 
and hence Eq. 3.7 is hyperbolic, parabolic or elliptic at 
a point according to whether the discriminant b*- 4ac is 
positive, zero, or negative. 

To determine the characteristic curves we return to the 
pair of first-order Eqs. 3.1 and suppose that the system is 
hyperbolic throughout the domain of interest. At every point 
there are two roots, (dy/dx)a and (dy/dx)B , to the quadratic 
Bos. 3.4. A curve, which at each of its points has the slope, 
(dy/dx)a , is said to be a a characteristic; A,curve, whose 
slope is everywhere, (dx/dy)8, is said to be a B characteristic. 
Thus there are two families of characteristic curves filling 


the domain as shown in Figure 3. 
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FIGURE 3 


ЈЕ we are considering a characteristic direction so that 

the determinant in Eq. 3.3 is zero, the right-hand column 
must be compatible with this if there are to be any solutions 
wall for the first derivatives; i.e., when the right- hand 
column is substituted for any of the columns on the left, 

the resulting determinant must also vanish. Replacing the 
fourth column on the left with the column on the right and 


setting the determinant equal to zero leads to 
dy - 
(A,B2-A7B])du т [A1 €2-^56,) pe (805-801). dv = 0. (3.8) 


Insert (dy/dx)a into Eq. 3.8; ‘it then becomes an ordinary 
differential equation of u and v along the characteristic. 


A similar equation can be obtained along a ß characteristic. 
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A method of attack has been outlined for solving hyper- 
bolic systems: first, locate the characteristic curves; 
second, integrate the ordinary differential Eq. 3.8 along 
the characteristic. This procedure is called the method 


ШЕ characteristics. 


B. FORMULATION OF THE EQUATIONS PROGRAMMED 

In Section II, D, Eq. 2.41 was developed as the partial 
differential equation for $(x,r) that needed to be solved 
subject to the boundary condition Eq. 2.42. If M> 1, the 
preceding equation is of the hyperbolic type and possesses 
real characteristics, which satisfy the ordinary differential 


equation 
(M? - l)dr? - dx? = 0 : (3.9) 


Introduce along the characteristics the arc-length ds^zM^dr^? 


с лас УМ ИН MCI n И! 3.10 
ЕТЕЙ опе Наг Х' = <= = == 7 т E = ( ) 
Note: The angle between the characteristics and the x-axis 
is Sina = 1 (3-11) 
M 


An arbitrary function F(x,r) has the following derivative 


along a characteristic 


ШЕ | i Е ‚ _ M*-1 y 
E c F. X Da r = Ue 





l 
M T M F. (3.02) 


Denoting the derivative along the left-running characteristic 


wath F the derivative along the right-running characteristic 


T: 
with Р, опе has 











/M2 1 1 M*-1 
F = -- = = Т 
1 ма = "и P. Mo С S F. (SUS) 
Solving for F,, F, results in | : 
FP. = E (p, - >.) M (FQ Fo) 
r 2 1 дА" К Е 1 2 (3.14) 


s Ум 


Са11 A the arc-length of the left-running characteristic, 


then one has for the second derivatives 


M 


ЕШ doe ; Fol = Tem (3.15) 
l 
Hence 
/M?-1 

Fig = (Fily * M _ (Fi), M 
and 

B u -M2 ae 

=> E. E E no 3 Prr | Sa 


Therefore the basic equation for $(x,r) was written 


2 
Emo, + = а 72 С. С +02) + кеме- 2, | Ф = 0 
Enc yo - Я 
(3.16) 
ФЕ 





-=È : m | 
TRA? 1KM 2E 


This equation was then written in finite difference form. 


Using the backward difference form of the finite dif- 


ference method as developed in Wang Ref. 3 one obtained 


оф | А | 
39) j 7s CEN +95] Applying this technique to 
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Eq. 3.16 where the points A and B were on the right-running 


characteristic (See Fig. 3) then 


me с=ш= шш AA ARA 
— -- 


AS 2rM ML 





ШЫН) - у Е Оо Eu (9, + 2 + («- n? 


FIGURE 4 


The values of Ф, 6, and Ф. in Eq. 3.17 were not specified at 
a point. When this equation was first written these values 
were taken at A, the previous point. It was found in 
comparing the results in Chapter V that there was a consider- 
able amount of error present. In investigating the source 

of these errors it was found that if the values for 9, $4 and 
Ф. Were averaged between points A and B that these errors 
disappeared. The value for ® at point B was not yet known 


but it was found by integrating along the characteristic. 


Using trapezoidal rule 


(B) = (A) + Fi Ф. (А) + гэ во] аа (3218) 


Substituting this into Eq. 3.17 where (B) arises, and using 


= 


the above averaging technique for %), %> and 9? results in 
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– 45 


тву AS? < I 
Ф. (В) 4r(B) ue 1 AS - Бл K ~- ВИ: 
(А) sais iKM AS? а | 
Ф» 4r(A)M * 2yM?^-] B9 - 3- (К - =a + (3.19) 


B n? AS n AS 
a (к трут.) 88 5 00А) (“- заву] ж 


аа тах Ту to solve for.95 at point B, the above procedure 
was used except now move along the left running character- 
istic and integrate 91 along the left running characteristic 


fo find Ф(В). 


AS iKM iKM 
pns 1 - dE(B)M " Sr AS = : мейе qe СӘРЕ SIs 48) т 


NS iKM 2 
n0 [som * Se 29 c P i^ -m»m]] - 





_AS ikM AS? n2 

$4 (C) [a em ая 2742-1 Se 4 (е ar) | + (3.20) 
ne AS n2 

|“ ке тома “е [к "em | 25 


Then there are two equations [3.19 and 3.20) in two 
unknowns, $, (B) and 9$4(B), which can be solved by using 


2 RE 
Cramer's rule. The following abbreviations are defined as 


- 2 
Crandall, S. H., Engineering Analysis, New York, 1956, 
Di 26. 
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a * AC AS 





Я — ÀS iKM 
EN ! “15114 - 2747-i ^? 

= E EE 
Маша” mem - mer ~ 
В = шара ate 1KM AS - 452 2 n? (3.21) 

Ar (B)M > M2=1 ok а E 
С = ф)(АЈАЈ+ ф,(А) Е E - A TT AS^ [ o2 n? 
1 n? l 2 
= 0 2 — —— 9 AN IS + + % pL 
- er ec SINT S2 n? 

ilo c» s | ancora” arma? 7 (e - wr) | 

2 (с) [xe- — n° = 2 ом 

> о EXC AS + 5 Ф(А) қ - тТ(Бум? 5 

In terms of these abbreviations the two equations are 

ВФ. (В) + А Ф, (В) = р (34225 
Then using Cramer's rule 

$, (B) = CA - BD/A* - B* 

DE pa а CD/A? - pg? (3228) 

94 was integrated to find 4galong the right running char- 
acteristic and as a check 0, was also integrated along the 


left running characteristic using the trapezoidal rule 


= 
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1 ú 1 | 
Ф (В) = Ф(А) + 5 | (А) + $ (B) | AS 
2 = de 
(В) - (C) + 5 D + s o) AS | (3,23) 
The difference |o!) - о? (в) | was a measure of the 


accuracy. If it becomes too large, the grid size must be 


made smaller. Ф! (Вв) and Ф? (В) were averaged together for 


the best results. 
Ф (В) = š |0: (в) + 3° (в) (3.25) 


The above procedure will suffice for a general field 
point which is not subject to a boundary condition once the 
procedure has been initiated. 

on the initial Mach line iv isSfknowni thot om O P ML LE 


дфу = 0, therefore from Eqs. 3.20 and 3.21 (See Eg ӘЙ 


= Ф 
E A, ТБ) 


$, (E) - DD TE (3.26) 


N М» 
E ж 


ES 2 
S ` 
ES ` 
S : 


` 
ETGURES>S 
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At the initial point on the initial Mach line the panel 


boundary condition also applied, i.e., 


UON ae E 
pH Е li іка] : T2 
puso from Eq. 3.14 


x 
DNE (a, - 9) 


remembering that C= 0 and combining the two equations 





2 7 
DM. xe аа E ь 
а М E: ixz] | ны 


On the panel from the above boundary condition (See 


Ша. 6) 


= ОИЕ; 
$) (L) Зак) | M E + іка | (3528) 


End from Eqs. 3.19 and 3.21 9, (1) + $5(L) was found on the 
right running characteristic and once again there exists a 


system of two equations in two unknowns. 


E = iz] (87 28) 
9х 


Ao, (L) + Bö, (L) = С 


óJ (L) e Фо(ђ) = = 


SIND 


Using Cramer's rule? 


2 [32 , i ))/ 
Фу (1) = сва Z + 3x0) A+B 

_ 2 faz ка) | / 
ф = A —— И A+B (3.30) 
E" | “м % 
с И | 
Шота. 
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FIGURE 6 


To write the program for a vibrating panel or the two- 
dimensional case, the radius of the cylinder r was set to 
infinity. Inspection of the basic equations programmed 
(Eq. 3.21) showed that r was always in the denominator, 
which causes these terms to equal zero. Therefore, terms 
containing r were neglected in the Fortran IV coins of 
these equations. For the vibrating shell or three dimensional 
case the full equations as developed were Fortran IV encoded. 

Also, it should be noted that the problem was formulated 
for an arbitrary wall deflection Z(x). However, in the 
numerical computations only the case of a sinsuoidal standing 


маме given by 


Z(x) = Zo Sin (mtx) (3.31) 
was considered. The axial wave-length l is hence equal to 
L/2m where m is the number of axial half-waves in a given 


cylinder length L. 
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"TEST CASES 
The first test case considered a two-dimensional flat 
panel. If the cylinder radius is set at a very large value 
in relation to the length, then one has essentially a flat 
panel in supersonic flow. | 
Te ne a 
R > >а 


When MFREQ is set at various values, NFREQ and K set at zero, 
one then has supersonic flow over a wavy wall as developed 


an Johns Ref. 1 . 


R >>| 


Taking the wavy wall equations for 6$, oF and Ф. from Johns 


Ref. 1 and normalizing results in 


ф = - ART Sin (x - vM*-1 r) 
АЛ Cos Е (x — /м2-1 ) | | 


юн 
|| 


е. = Tic Cos E - /M*-1 9] (32.51) 
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Ву 
it 
to 
to 


to 


using a large value of R in the three-dimensional case, 
and the two-dimensional case results were then compared 
the wavy wall results. These results were found to agree 
the fifth significant figure. This case was also used I 


remove the majority of minor programming errors. 


To remove any other hidden programming errors it was 


decided to set each variable equal to a finite value and to 


step through the first three characteristic lines by hand. 


The values were then compared to those computed by the 


computer. In this manner the rest of the programming bugs 


were found and removed. 
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IV. COMPUTER PROGRAM ORGANIZATION AND OPERATION 


A. GENERAL ORGANIZATION 

The numerical methods discussed in the foregoing 
sections have been coded as a FORTRAN IV computer program 
for the IBM 360 system. 

The program consists of a main control section and four 
subsections or subroutines. An algorithm showing the 
program information flow and the interaction between the 
control section and the subroutines was listed in Figure 7. 
The control section reads in the input data, computes the 
characteristic grid, computes the initial condition, and 
computes several constants common to the subroutines. Then 
subroutine COMPRX is called to compute the а 
of the next grid point to be calculated, see Figure 8 for 
a diagram of the characteristic grid used. Two "logical if" 
statements then determine if the point in question is 
located on either the initial Mach line, the panel, or at 
a general field point. Control is then passed to that 
subroutine to compute the velocity potential and its 
derivatives along the characteristics at that point. The 
coefficient of pressure is computed along the panel in 
subroutine PANPT and after the program has stepped down to 
Pie last point on the last characteristic, a listing of the 


coefficient of pressure at each panel point is printed out. 
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COMPUTE 
READ CONSTANTS, 


SET 
INPUT З 
COUNTERS 





WRITE 
INITIAL 
VALUES 






SET LINE 
CALL COUNTER 
COMPXR = 2 


CALCULATE 
ICS 









COMPUTE 
CONSTANTS 


ONSTANTS 





Yes 


CALL 
PANPI 








POINT ON PANEL 


NO 


CALL 
GENFPT 


Yes 







WRITE VALUES END 














OF OF FINAL RESET 
Мо. Ф CHARACTERISTIC FLAGS & 
| 1 2 LINE COUNTERS 





Yes 


Yes 
COEFFICIENTS 


FIGURE 7. 
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B. MAIN CONTROL PROGRAM OPERATION 

The main program dimensions the storage requirements 
in common block form; then it reads three input data cards. 
Кит 5- сака has the date of the run in Paro ee 
second card reads in the fineness of the grid, number of 
this run, circumferential mode number, and the axial mode 
number in 415 format. The third card reads Mach number, 
the reduced frequency number, and the cylinder radius in 
BEFO.4 format. Input is then written to verify input data: 

The main program then computes the Mach angle, step 
Sizes, and zeros the counters. The values for the velocity 
potential and its derivatives and the coordinates of each 
point are stored aS a nx 2 matrix. As the program steps 
from characteristic: line to characteristic line, only those 
values for the line being computed and the line just 
computed are in storage. The initial values at the first 
point on the initial characteristic is computed using Eq. 
3.27 and written out to aid in the program check-out. The 
line counter, "Kount", is then set at two and subroutine 
COMPXR is called to compute the values of x and r at the 
next point. The main program next computes various quantities 
dependent on x and r which are needed in the other subroutines. 
The values for x, r, and dependent quantities are printed out 
as a program debug aid. Two "logical if" statements next 
determine if the point under consideration is on the initial 
Mach line, at a general field point, or on the panel. The 
appropriate subroutine is then called to solve for the 7 


velocity potential and its derivatives along the 
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Shuedeteristics at that posnt. The next оаза са не 
statement determines if the program must switch lines 
Step along the present line. Two оге Иос еса а 
statements are used to determine if this solution is 
completed and if there is more data to be read and more 
solutions to be computed. If "numrun" is set at two, 

the program will look for more data; if "numrun" equals 
one, the program will terminate. Before it stops or 
reads more data the coefficient of pressure along the 
panel will be printed out for each panel point and an in- 
line graph, using the IBM scientific library program POLTP, 
will plot the real and imaginary parts of the coefficient 


Of pressure. 


e SUBROUTINE COMPXR OPERATION 

Information is passed between the calling program and 
the subroutines by means of the common storage. The first 
step in COMPXR is to determine if the point to be computed 
is on the initial Mach line or not. The coordinates of the 
last point computed are used and the step distance in the 
x and r directions either added to or subtracted from them 
to find the new x and r coordinates. The addition or sub- 
traction depends upon whether we step up to a new characters 


istic line or down the present characteristic line. 
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D. SUBROUTINE MACHLN OPERATION 

In subroutine MACHLN the boundary conditions of the 
equations developed in the foregoing section are enforced 
such that the velocity potential and its derivatives along 
F right running characteristic are set to zero. есца сане 
3.26 were coded in FORTRAN to solve for 95 at the point in 
question. “he value of $, and the хатари исса собе 
Bor? are printed out for the first 12 characteristic е6 
to aid in the program checkout and then control is returned 


to the main program. 


в, OPERATION OF SUBROUTINE PANPT 

In subroutine PANPT the boundary condition along the 
panel is enforced. The velocity potential and its derivatives 
along the characteristics are computed at the Set in 
question, using Eg. 3.29. The real and imaginary parts of 
the coefficient of pressure are stored for future read out. 


Control is then returned to the main program. 


F. OPERATION OF SUBROUTINE GENFPT 

GENFPT computes the velocity potential and its first 
derivatives along the right and left running characteristics. 
It draws the necessary information from common storage and 
Mees the Eqs. 3.21, 3.22, and 3.23, which were аечетореа п 
EHe foregoing chapter. To aid in program check out the 
variables A, Ај, A5, B, C, and D are printed out for the 
first 12 characteristic lines. Control is then returned to 


the main program. 
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V. RESULTS AND DISCUSSTON 


The case of a vibrating two-dimensional panel was 
chosen as the first test of the accuracy of the character- 
istics method. These results were compared to the results 
obtained by evaluating the method-of-singularities integrals 
developed in Ref. 4 for flutter of a two-dimensional panel 
with one surface exposed to supersonic flow. These integrals 
were evaluated by applying a 12-point Gaussian Quadrature. 
The results of these evaluations were available in unpublished 
notes by M. F. Platzer. 

The comparisons for M - v2 at a reduced frequency of 0.5 
and an axial mode number of one are shown in Figures 8-11. 
Figures 8 and 9 ea of the real and imaginary parts of 
the coefficient of pressure obtained by the method of singu- 
larities. Figures 10 and 11 are plots of the real and 
imaginary parts of the Cp for the characteristics method. 
When these tests were first run it was found that the real 
parts of the Cp plots were identical, as shown in Figures 8 
and 10. The imaginary parts of Cp, shown in Figures 9 and ll, 
did not agree. In rechecking the formulation of the equations 
it was found that the values of Ф, 2, and 5 had not been 
averaged as explained in Chapter II, Section B, where Eqs. 
3.19 and 3.20 were formulated. When these average values 
were incorporated into the equations, the computer program 
rewritten and the program rerun, then the two plots were 


found to be identical for both the real and imaginary parts 


of the Cp. 
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Figures 12-14 are plots for a reduced frequency of two, 
an axial mode of four, and the same Mach number. These 
comparisons of the real and imaginary parts of the Cp were 
also found to be identical. 

In attempting to verify the results of the vibrating 
shell program, Ref. 5 was used. Reference 5 includes a 
study of supersonic potential flow past an infinitely long 
Vibrating cylinder. It was not known what effect the dif- 
ference in an infinitely long cylinder vice a finite length 
cylinder would have on the outcome, but these were the only 
published results found that could be used for a comparison 
of the vibrating shell program. 

Two parameters were used to compare the results of the 
two studies: the maximum amplitude, Ај, of the sinsuoidal 
Cp wave, and the phase angle, ў], of this wave. Equation 14 
in Ref. 5 is used to equate Cp to the value of А, listed in 


Figure 2a of Ref. 5. 


= 2 U? mıR mTx 
Bx,R,9) - Er -= T Cos ng Ay Cos E ar 01 





р is defined as p = p ~ Po where 
p 7 po 
= Loy? 7 Sc 
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and the terms, Z , L, and Cos(x) were all taken as equal 


о” 
ШЕ опе then 
С = Do == 2 
р ^ pU M BUE 
M 
Ai * Cp om т. (5.2) 


Three different values of n were computed and the coefficients 
Of pressure were listed in Figure 17 as a test case for this 
comparison. For the case n = 0, (See Figure 2a of Ref. 5) 

we find that for the ratio of axial wave length/cylinder 
Bagius 2/R = 3.33, A, equals 1.05. Using Eq. 5.2 to determine 


E For a maximum Cp of 5.63 in Figure 17, Ај for this study 


would equal 1.044. Using this same procedure for п = 4 and 8 
n = 4 A, (ref 5) = 1.30 A, = 1.445 
п = 8 A, (ref 5) = 1.80 Ai = 1.85. 


TO compare the phase angle Vi» Figure 2b in Ref. 5 was 
entered for 2/R = 3.33 and the three values of n taken 
previously. This shift, ф ү, was calculated from results 


(Fig. 17) by determining the zero crossing point. 


n = 0 59.8 
n = 4 БИЯ 
n= 8 48.9 


Converting these zero crossing points to angular shifts, 


me was found that 
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n= 0 V1 = 4.8 расе ў 
п = 4 V1 = 13.5 Ui (ref 5) = 10 
n= 8 МА = 55.4 Ф, (кеЁ 5) = 60 


which compare favorably with Anderson's values listed in 


mae third column. 
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SHILIYVINONIS JO AOHLUN 
TANVd јула IVNOISNSHIG-OML DNILVHGIA V НОЯ dO 40 LUVA AXUVNIVOVNI 


Ol Yunald 


g' 9” ЛИХ b° >: 


S'O =ġ' W'A =N v 
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SOILSIUHLOVUVHD JO AOHLSN 
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SYILIUVININIS dO GOHLYN 


TANVd LVId IVNOISNENIG-OML DNIGVEGIA WV НОЯ do JO LHVd 'IVWH 


ET JENDIA 
1/х 


do Jo LUYd IVAU 
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SHILZLIEVINONIS JO GOHISN 
яуа IVIA TVNOISNANIA-OMI ONILVEEIA У НОЯ ао мо Luvd AYVNIVOVWI 
TL ЗЯПЭІЧ 


ay 
СІ 8' 9 t° аў 


do JO інуа XEVNIDVHI 
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01 8: 2 > e о 


do до дама ТУН 
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ОТ FUNDI 


T/R 
О] 8” 9 v с О 


1 q =u zur 


N 
| 


do JO тиха AEBVNIOVWI 
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COEFFICIENT OF PRESSURE COMPUTER OUTPUT LISTING 


FOR CASE OF VIBRATING SHELL 


паэ Э, м = 3, ВО, Е 0.20 and n = 0, 1, 3, 


where L is divided into 120 grid points. 











Grid 

Point n= 0 n= 4 n= 8 
lh. 5.6198 5.6198 5.6198 
2. 5.5677 5.5538 5.5642 
3. 5.4813 5.4266 5.4676 
4. 5.3613 5.2396 5.3307 
5. 5.2083 4.9950 5.1546 
6. 5.0233 4.6954 4.9404 
T 4.8075 4.3440 4.6897 
8. 4.5621 3.9444 4.4041 
9. 4.2887 3.5007 4.0857 
10 3.9890 3.017? 3.7366 
BL. 3.6647 2.4988 3.3591 
12 3.3180 1.9504 2.9559 
13 2.9509 1.3775 2.5295 
14. 2.5657 0.7855 2.0828 
15. 2.1648 0.1801 1.6188 
16. 157507 =0. 4329 1.1407 
17. 1.3258 -1.0477 0.6515 
183. 0.8929 -1.6586 0.1546 
19. 0.4546 -2.2597 -0.3467 
20. 0.0136 -2.8453 -0.8492 
Dl. -0.4273 -3.4100 -1.3495 
22. -0.8656 -3.9435 -1.8442 
po. -1.2985 -4.4557 -2.3301 
24. -1.7232 -4.9267 -2.8040 
25. -2.1373 -5.3572 -3.2626 
26. -2.5380 -5.7431 -3.7028 
27. -2.9230 -6.0806 =4 loo 
28. -3.2899 -6.3665 -4.5167 
29. -3.6365 -6.5980 -4.8848 
30. -3.9605 -6.7726 -5.2237 
31. -4.2600 -6.8886 -5.5309 
32. -4.5331 -6.9446 -5.8044 
33. -4.7782 -6.9398 -6.0423 
34. -4.9938 -6.8738 -6.2428 
35. -5.1785 -6.7468 -6.4046 
36. -5.3312 -6.5595 -6.5263 
31. -5.4509 -6.3132 -6.6071 
38. -5.5369 -6.0095 -6.6462 
39. -5.5887 -5.6507 -6.6432 
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п = 4 





6.6619 
6.3383 
OD 
5.5841 
2 DIG 
4.7026 
4.2205 
3.7146 
3.1878 
2.6436 
2.0848 


ISSO 


0.9376 

0.3559 
-0.2264 
-0.8061 
-1.3795 
-1.9434 
-2.4942 
-3.0288 
-3.5439 
-4.0364 
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-4.9423 
-5.3503 
-5.7250 
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VI. SUMMARY AND CONCLUSIONS 


A linearized method of characteristics procedure was 
developed to compute supersonic flow past vibrating panels 
and shells. Good agreement was obtained when comparing 
with previous results by Nelson and Cunningham [Ref. 4] 
and Anderson (Ref. БІ? who used quite different approaches 
to this problem. 

However, more work needs to be done to further verify 
the cylindrical shell computations. In particular, the 
computations of the generalized aerodynamic forces presented 
by Dowell and Widnall (Ref. 6 | should be verified. For this 
purpose the limiting case of vanishing shell radius (slender 
body theory) must be studied. This requires a ren on ата 
of the flow boundary condition. 

Also, it should be noted that the characteristics method 
allows the incorporation of quite arbitrary vibration modes 
(not restricted to sinsuoidal modes) and hence should prove 
to bea quite versatile tool for aeroelastic analyses in the 


supersonic flow regime. 
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